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We describe early success in the evolution of binary black hole spacetimes with a numerical code 
based on a generalization of harmonic coordinates. Indications are that with sufficient resolution 
this scheme is capable of evolving binary systems for enough time to extract information about the 
orbit, merger and gravitational waves emitted during the event. As an example we show results 
from the evolution of a binary composed of two equal mass, non-spinning black holes, through a 
single plunge-orbit, merger and ring down. The resultant black hole is estimated to be a Kerr black 
hole with angular momentum parameter a « 0.70. At present, lack of resolution far from the binary 
prevents an accurate estimate of the energy emitted, though a rough calculation suggests on the 
order of 5% of the initial rest mass of the system is radiated as gravitational waves during the final 
orbit and ringdown. 



/. Introduction: One of the more pressing, unsolved 
problems in general relativity today is to understand 
the structure of spacetime describing the evolution and 
merger of binary black hole systems. Binary black holes 
are thought to exist in the universe, and the gravitational 
waves emitted during a merger event are expected to be 
one of the most promising sources for detection by grav- 
itational wave observatories (LIGO, VIRGO, TAMA, 
GEO 600, etc.). Detection of such an event would be 
an unprecedented test of general relativity in the strong- 
field regime, and could shed light on many issues related 
to the formation and evolution of black holes and their 
environments within the universe. Given the design- 
goal sensitivities of current gravitational wave detectors, 
matched filtering may be essential to detect the waves 
from a merger, and extract information about the as- 
trophysical source. During the early stages of a merger, 
and the later stages of the ringdown, perturbative ana- 
lytic methods should give a good approximation to the 
waveformj^ 1^; however, during the last several orbits, 
plunge, and early stages of the ring down, it is thought a 
numerical solution of the full problem will be needed to 
provide an accurate waveform. 

Smarrj^ pioneered the numerical study of binary black 
hole spacetimes in the mid-70's, where he considered the 
head-on collision process in axisymmetry. The full 3D 
problem has, for many reasons, proven to be a more 
challenging undertaking, and only recently has progress 
been made in the ability of numerical codes to evolve 
binary systems 0, Q IE IE Hi- However, until now no 
code has been able to simulate a non-axisymmetric col- 
lision through coalescence and ringdown. The purpose 
of this letter is to report on a recently introduce numeri- 
cal method based on generalized harmonic coordinates 
that can evolve a binary black hole during these crucial 
stages of a merger. At a given resolution the code will 
not run "forever" , though convergence tests suggest that 
with sufficient resolution the code can evolve the system 
for as long as needed to extract the desired physics from 
the problem. As an example we describe an evolution 
that completes approximately one orbit before coales- 
cence, and runs for long enough afterwards to extract 



a waveform at large distances from the black hole. 

The code has several features of note, some or all of 
which may be responsible for its stability properties: (1) 
a formulation of the field equations based on harmonic 
coordinates as first suggested in (2) a discretization 
scheme where the only evolved quantities are the covari- 
ant metric elements, harmonic source and matter func- 
tions, thus minimizing the number of constraint equa- 
tions that need to be solved (which is similar to the dis- 
cretization scheme used in 11]), (3) the use of a com- 
pactified coordinate system where the outer boundaries 
of the grid are at spatial infinity, hence the physically 
correct boundary conditions can be placed there, (4) the 
use of adaptive mesh refinement to adequately resolve 
the relevant length scales in the problem, (5) dynam- 
ical excision that tracks the motion of the black holes 
through the grid, (6) addition of numerical dissipation 
to control high-frequency instabilities, (7) a time slicing 
that slows down the "collapse" of the lapse that would 
otherwise occur in pure harmonic time slicing, and (8) 
the addition of "constraint-damping" terms to the field 
equations 01 . This final element was not present in 
the version of the code discussed in and though these 
terms seem to have little effect when black holes are not 
present in the numerical domain, they have a significant 
effect on how long a simulation with black holes can run 
with reasonable accuracy, at a given resolution. 

An outline of the rest of the paper is as follows. In 
Sec. II we give a brief overview of the numerical method, 
focusing on details not present in 9J. Sec. Ill gives re- 
sults from the simulation of one such orbital configura- 
tion. We conclude in Sec. V with a summary of future 
work. More details, including convergence tests, the ef- 
fect of constraint damping and a thorough description of 
the initial data calculation, will be presented elsewhere. 

II. Overview of the numerical method: We briefiy 
summarize the formulation of the field equations, gauge 
conditions and initial data used here, emphasizing de- 
tails that are not contained in We discretize the 
Einstein field equations expressed in the following form 
(using units where Newton's constant and the speed of 
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light are equal to 1): 

g^^gaP.-iS + g'^^ ,l39aS,'y + g'^^ Mgi3S,f + '^H(^a,l3) 

-K{naCi3 +nf3Ca - ga0n''Cj) ■ (1) 

are (arbitrary) source functions encoding the gauge- 
freedom of the solution, F*^ are the Christoffel symbols, 
Tap is the matter stress tensor with trace T, k is a pos- 
itive constant multiplying the new constraint damping 
terms followingH^, n'' = l/a{d/dtY - P'/aid/dx'Y is 
the unit hypersurface normal vector with lapse function 
a and shift vector {x^ = i, = [x-'^, x^, x^] = [x, y, z]), 
and are the constraints: 

Cf, = Hf,^gf,,Dxr (2) 

We use the following to evolve the source functions: 

DHt = -6^^ + 6i^t..n^ H,^0 (3) 
a'' 

where and ij are positive constants. Note that Q is 
not the usual definition of spatial harmonic gauge, which 
is defined in terms of contravaricnt components H^. 

We use scalar field gravitational collapse to prepare ini- 
tial data that will evolve towards a binary black hole sys- 
tem. Specifically, at t = we have two Lorentz boosted 
scalar field profiles, and choose initial amplitude, sepa- 
ration and boost parameters to approximate the kind of 
orbit that the black holes (which form as the scalar field 
collapses) will have. The procedure used to calculate the 
initial geometry is based on standard techniauesfl^. and 
is a straight forward extension of the method described 
in|^ to include non-time-symmetric initial data. The ini- 
tial spatial metric and its first time derivative is confor- 
mally flat, and we specify a slice that is maximal and 
harmonic. The Hamiltonian constraint is used to solve 
for the conformal factor. The maximal conditions K ^ 
and dtK = {K is the trace of the extrinsic curvature) 
give the initial time derivative of the conformal factor 
and an elliptic equation for the lapse respectively. The 
momentum constraints are used to solve for the initial 
values of the shift vectors, and the harmonic conditions 
— are used to specify the initial first time deriva- 
tives of the lapse and shift. 

Results: In this section we describe results from the 
evolution of one example of a scalar field constructed bi- 
nary system. The present code requires significant com- 
putational resources to evolve binary spacetimes|l8||. and 
thus to study the orbital, plunge, and ringdown phases 
of a binary system in a reasonable amount of simulation 
time we chose initial data parameters such that the black 
holes would merge within roughly one orbit — see Fig. ^ 
and Table ^ The following evolution parameters in 
and (jSJl were chosen: k « 1.25/Mo, « 19/Mo, ^2 ~ 
2.5/Mo, ?7 = 5 (these parameters did not need to be 
fine tuned), where Mq is the mass of one black in the 
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FIG. 1: A depiction of the orbit for the simulation described 
in the text (see also Table||J. The figure shows the coordinate 
position of the center of one apparent horizon relative to the 
other, in the orbital plane 2 = 0. The units have been scaled 
to the mass Mq of a single black hole, and curves are shown 
from simulations with three different resolutions. Overlaid on 
this figure are reference ellipses of eccentricity 0, 0.1 and 0.2, 
suggesting that if one were to attribute an initial eccentricity 
to the orbit it could be in the range — 0.2. 



binary. This system was evolved using three different 
grid hierarchies, which we label as "low", "medium" and 
"high" resolution. The low resolution simulation has a 
base grid of 32-^, with up to 7 additional levels of 2 : 1 
refinement (giving a resolution in the vicinity of the black 
holes of ~ Mq/IO). For computational efficiency we only 
allowed regridding of level 6 and higher (at the expense 
of not being able to accurately track out-going waves). 
For the medium resolution simulation, we have one ad- 
ditional level of refinement during the inspiral and early 
phases of the merger, though have the same resolution 
over the coarser levels and at late times; thus we are 
able to resolve the initial orbital dynamics more accu- 
rately with the medium compared to low resolution run, 
though both have roughly the same accuracy in the wave 
zone. The high resolution simulation has up to 10 lev- 
els of refinement during the inspiral and early ringdown 
phase, 9 levels subsequently, and the grid structure of the 
lower levels is altered so that there is effectively twice the 
resolution in the wave zone. The reason for this set of hi- 
erarchies is again for computational efficiency: doubling 
(quadrupling) the resolution throughout the low resolu- 
tion hierarchy would have required 16 (256) times the 
computer time, which in particular for the higher resolu- 
tion simulation is impractical to do at this stage. 

Fig. 121 shows the horizon masses and final horizon an- 
gular momentum as a function of time. The ADM mass 
of the space time suggests that approximately 15% of the 
total scalar field energy does not collapse into black holes. 
The remnant scalar field leaves the vicinity of the orbit 
quite rapidly (in t « 30Mo, which is on the order of the 
light crossing time of the orbit). Black hole masses are 
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low res. 


med. res. 


high res. 


ADM Mass 


2.36Mo 


2.39Mo 


2.39M0 


initial BH masses 


0.97Afo 


0.99Mo 


Mo 


orbital eccentricity 


0-0.2 


0-0.2 


0-0.2 


proper separation 


16.5Mo 


I6.6M0 


I6.6M0 


angular velocity xMq 


0.023 


0.023 


0.023 


final BH mass 


1.77Afo 


1.85Mo 


I.9OM0 


BH spin parameter 


0.74 


0.73 


0.70 



TABLE I: Some properties of the simulated equal mass binary 
system described in the text. Where relevant the units have 
been scaled to the mass Mq of one of the initial black holes, 
measured from the higher resolution simulation at a time after 
the majority of scalar field accretion has occurred. The final 
black hole mass and spin where estimated from data as shown 
in Fig. 121 a little while after the black hole formed, though not 
so long after as to be affected by the "drift" from numerical 
error. The initial proper separation was measured at t — 
lOMo, and is the proper length of the piece of a coordinate line 
outside the apparent horizons that connects their coordinate 
centers. The black holes initially have zero spin. 
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FIG. 2: The plot to the left shows the net black hole mass 
of the spacetime in units of the mass Mo of a single initial 
black hole, calculated from apparent horizon (AH) properties 
(using Q with the dynamical horizon estimate for J), and 
from simulations with three different resolutions. The initial 
sharp increase in mass is due to scalar field accretion, the 
small "wiggle" at around 20A/o appears to be a gauge effect, 
and the "jaggedness" around the time of the merger is due 
to robustness problems in the AH finder that manifest when 
the AH shapes are highly distorted. To the right the Kerr 
parameter a of the final black hole is shown (for clarity we 
only plot the results from a single simulation), calculated us- 
ing the ratio Cr of polar to equatorial proper circumference 
of the AH and applying (O, and using the dynamical hori- 
zon framework (curve labeled DH). The loss of mass (and 
similarly increase in a) with time after the merger is due to 
accumulating numerical error. 



estimated from the horizon area A and angular momen- 
tum J, and applying the Smarr formula: 



M 



^Ml + jy{AMl), (4) 



The horizon angular momentum of the final black hole is 
calculated using two methods (which do give zero angu- 
lar momentum when applied to the initial black holes. 



as expected). First, by using the dynamical horizon 
frameworkjlSj. though assuming that the rotation axis 
of the black hole is orthogonal to the z = orbital plane, 
and that each closed orbit of the azimuthal vector field 
(which at late times should become a Killing vector) lies 
in a z = constant surface of the simulation. Due to the 
symmetry of the initial data, these assumptions are prob- 
ably valid, though this will eventually need to be con- 
firmed. The second method, following 16], is to measure 
the ratio of the polar to equatorial proper radius of the 
horizon, and use the formula that closely approximates 
the function that is valid for Kerr black holes: 



1 - (2.55a - 1-55)^ 



(5) 



As seen in Fig. |21 the initial ringing of the black hole 
is quite apparent in the estimate using Cr- Remarkably, 
the dynamical horizon estimate for a and average value 
obtained using Cr agree quite closely, even shortly after 
the merger when one might have expected the black hole 
to still be too far from its stationary state to have either 
method be applicable. 

To estimate the gravitational waves emitted by the bi- 
nary we use the Newman- Penrose scalar ^'4, with the null 
tetrad constructed from the unit timelike normal n'^, a 
radial unit spacelike vector normal to r = constant coor- 
dinate spheres, and two additional unit spacelike vectors 
orthogonal to the radial vector^^. Far from the source, 
the real and imaginary components of ^4 are propor- 
tional to the second time derivatives of the two polar- 
izations of the emitted gravitational wave. Fig. |21 shows 
an example of the real part of ^'4. Most of the early, 
short wavelength burst of waves can be correlated with 
the passage of the remnant scalar field that did not fall 
into the black holes (the "noisy" nature of this piece of 
the waveform is in part due to numerical error). This 
unwanted radiations leaves the domain quite early on, 
and so does not significantly affect the subsequent merger 
waves. Roughly the first long wavelength oscillation in 
the plot can be associated with orbital motion, and sub- 
sequent waves with the ringdown of the final black hole. 

To estimate the total energy E emitted in gravitational 
waves, we use the following formula |0| 



dE [ ,^ f\ , 

— = — / pdil, P= '^4dt ■ 



dt 



An 



*4di, (6) 



where ^4 is the complex conjugate of ^'4, and the surface 
integrated over in (jnj is a sphere of constant coordinate 
radius R (in uncompactified coordinates). This method 
of calculating the energy is quite susceptible to numeri- 
cal error, as we are summing a positive definite quantity 
over all time to give a change of energy with respect to 
time; thus numerical error in '54 will tend to inflate the 
answer. To reduce some of this error, we filter out the 
high spherical harmonic components {> ( — \m\ — 6) of 
'E'4 before applying ©. Note that the smaller integra- 
tion radii (as shown in Fig. ^ are not very far from the 
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FIG. 3: A sample of the gravitational waves emitted during 
the merger, as estimated by the Newman-Penrose scalar 4'4 
(from the medium resolution simulation) . Here, the real com- 
ponent of \['4 multiplied by the coordinate distance r from the 
center of the grid is shown at a fixed angular location, though 
several distances r. The waveform has also been shifted in 
time by amounts shown in the plot, so that the oscillations 
overlap. If the waves are measured far enough from the cen- 
tral black hole then the amplitudes should match, and they 
should be shifted by the light travel time between the loca- 
tions (i.e. by 25Afo in this example) . That we need to shift the 
waveforms by more than this suggests the extraction points 
are still too close to the black hole; the decrease in amplitude 
is primarily due to numerical error as the wave moves into 
regions of the grid with relatively low resolution. 

binary system, and so possibly in a region where lO is 
not strictly valid. However, the larger integration radii 
are in regions of the grid that do not have very good 
resolution (due both to the mesh refinement structure 
and the spatially compactified coordinate domain), and 
so numerical error (mostly dissipation) tends to reduce 



the amplitude of the waves with distance from the source. 
With all these caveats in mind, the numbers we obtain 
from © are 4.7%, 3.2%, 2.7%, 2.3% at integration radii 
of 25A/o, 50A'/o, 75Mo and lOOAfo respectively (from the 
high resolution simulation ,.2Qj), and where the percent- 
age is relative to 2Mq. Another estimate of the radiated 
energy can be obtained by taking the difference between 
the final and initial horizon masses (Table P) — this sug- 
gests around 5% (high resolution case). 

V. Conclusion: In this letter we have described a nu- 
merical method based on generalized harmonic coordi- 
nates that can stably evolve (at least a class of) bi- 
nary black hole spacetimes. As an example, we pre- 
sented an evolution of a binary system composed of non- 
spinning black holes of equal mass Afp, with an initial 
proper separation and orbital angular velocity of approx- 
imately IQ.QMq and 0.023/A/o respectively. The binary 
merged within approximately 1 orbit, leaving behind a 
blackhole of mass Mf « I.QMq and angular momentum 
J « 0.70M|. a calculation of the energy emitted in 
gravitational waves indicates that roughly 5% of the ini- 
tial mass (defined as 2Mo) is radiated . Future work 
includes improving the accuracy of simulation (in par- 
ticular the gravitational waves), exploring a larger class 
of initial conditions (binaries that are further separated, 
have different initial masses, non-zero spins, etc.), and 
attempting to extract more geometric information about 
the nature of the merger event from the simulations. 
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